Ab initio studies of opto-electronic excitations in VO2 
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We study the optical response of VO2 in the Ml insulating phase using methods based on density 
functional theory in its most recent developments. We start from a hybrid functional approach 
which may be a good starting point to carry out many-body perturbation theory since, as we show, 
it gives a qualitatively and to some degree quantitatively correct description of the insulating phase. 
In order to calculate the dielectric function, first, using hybrid density functional wavefunctions we 
have added GW corrections on top of the hybrid density functional calculations and we solved the 
so-called Bethe Salpeter equation to include effects of correlations of the electron-hole pair created 
upon photon absorption. We find that the effects of electron correlations are very important and 
they show up as a strong contribution of the electron-hole interaction in calculating the effects of the 
Bethe-Salpeter equation in the electron/hole pair excitation on top of the interacting group state. 



I. INTRODUCTION 

A plethora of fascinating new phenomena has been re- 
cently reported on oxide hetero-structures of transition 
metal oxides (TMO) and devices^— For example, an in- 
terface between two insulators behaves as a metal^ which 
becomes superconducting at sufficiently low tempera- 
tures, while an interface between two antiferromagnets 
becomes ferromagnetic^ These new structures not only 
create a playground for unexpected physical phenomena 
to be observed, but, in addition, they open up the pos- 
sibility for new applications based on radically different 
foundations. Probing such interesting new phenomena 
is quite difficult due to impurity effects, lattice imper- 
fections, and other materials-growth limitations. Very 
recent progress in material synthesis however, is 
starting to make it possible to carefully investigate these 
fascinating systems. 

The complex, unusual, and as yet not fully discovered 
or understood behavior of these strongly correlated ma- 
terials, such as transition metal oxides (TMO), can be 
manipulated in a variety of fundamentally new applica- 
tions. In particular, we would like to mention one such 
possibility which has been recently proposed 1 ^ by one 
of the authors of the present paper. This is based on 
the fact that these strongly correlated localized electrons 
form an electronic system which can be near the metal to 
Mott-insulator transition^ Photovoltaic devices based 
on carefully chosen doped Mott insulators can produce 
a significant photovoltaic effect. More importantly, it 
was found 1 ^ that if the Mott insulator is appropriately 
chosen, the photovoltaic effect can lead to solar cells of 
high efficiency, where a single solar photon can produce 
multiple electron/hole pairs, an effect known as impact 
ionization, in a time-scale shorter than the time charac- 



terizing other relaxation processes. Increase of solar cell 
efficiency due to this process has been proposed previ- 
ously for band-gap semiconductors; however, the effect 
is not significant there, because the characteristic time- 
scale for impact ionization is comparable to the time- 
scale for electron-phonon relaxation. The reason that 
this is not expected to hold for a Mott insulator is that 
the large Coulomb repulsion present in a Mott insulator 
leads to a large enhancement of the impact ionization 
rate. 

It is very important to have a computational ab initio 
scheme which is reliable to evaluate opto-electronic prop- 
erties such as excitations, gaps, absorption and, ideally to 
be able to compute the induced photo-current for a given 
bulk or interface structure. Furthermore, this scheme 
should be reliable for transition metal oxides where the 
electronic Coulomb correlations are expected to play an 
important role. Assuming that such a scheme exists, it 
would be very valuable in a variety of applications in real 
materials, including in the case of photovoltaic applica- 
tions proposed in Ref. Hol . This tool would meaningfully 
aid experimental efforts towards finding the most suit- 
able TMO based material with suitable size gap, band 
structure, size of transition matrix elements including se- 
lection rules, etc. 

In this paper we choose to study the optical properties 
of the insulating Ml phase of VO2. Vanadium dioxide 
is regarded as the prototypical example of a strongly 
correlated TMO material with a long history starting 
with the original suggestion by Mott himself^ and his 
collaborator Zylbersztejn nearly four decades ago. Fur- 
thermore, this material has been the playground for al- 
most all many-body techniques available to tackle strong 
Coulomb correlations. 

VO2 undergoes a structural distortion below approxi- 
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FIG. 1. (Color online) Unit cell of metallic and Ml phase of 
VO2 single crystal. The symmetry axes of these crystals are 
shown. The rutile c-axis is applied as a "parallel" direction for 
the distorted Ml phase in the context which is conventionally 
called a-axis in Ml phase. The experimental lattice constants 
are also depicted. 



mately T$ <~ 340-K" at ambient pressur e) 13 ' 14 . Above T$ 
the structure of VO2 is rutile- type and it is metallic and 
below Ts VO2 is an insulator. The insulating low tem- 
perature phase is monoclinic and is called the Ml phase. 
We would like to focus our studies on this phase because 
it is an insulator that exists at room temperature which 
could be utilized for photovoltaic applications. 

Below Ts and at room temperature, the vanadium (V) 
atoms dimerize and the V-V pair tilt around the rutile 
c-axis, doubling the unit cell as illustrated in Fig. [TJ This 
dimerization causes a zigzag-like antiferroelectric tilt of 
the V atoms perpendicular to this axis. The dimerization 
affects mainly the V 3d ti g states, which in the metallic 
phase fall near the Fermi energy and have similar band 
occupations ! 14 ' 15 In the dimerized state, these states hy- 
bridize differently to form dy and 7r* states. The bonding 
d\\ states fall below the ir* states, and they become fully 
occupied by the single d electron, leading to the insu- 
lating behavior due to this Peierls-type distortion. Our 
calculations support this picture and we illustrate what 
actually happens as the system undergoes the phase tran- 
sition from the metallic to the insulating Ml phase in 
Sec. MI O of this paper. 

Previous ah initio calculations of optical properties, 
such as of the real and imaginary parts of the zero mo- 
mentum dielectric function, show vast disagreement be- 
tween the result of the calculation and the experimental 
results. For example, dielectric results from the local den- 
sity approximation (LDA) within density functional the- 
ory (DFT^iii must be shifted entirely by hand to begin 
to agree with experiment^ Results from a generalized 
gradient approximation of DFT (GGA) together with a 
Hubbard-correction^ GGA+C/ look strikingly different 



apart from matching the optical gap2£ (which is effec- 
tively set by hand by tuning U). In addition, the elec- 
tronic structure is predicted to be metallic rather than 
insulating using both LDA and GGAi 21 ' 22 

Methods that go beyond ground-state DFT are now 
well established. However one should remember that 
a better starting point is absolutely necessary for d- 
clcctrons. 

With LDA+J7 one gets antiferromagnetic insulating 
ground states for both the Ml and the metallic phased 3 - 
In addition, it treats differently the various orbitals and 
we do not always a priori know which orbitals need this 
special treatment. Furthermore the exact value of the 
parameter U is not known and it is phenomcnologically 
determined. 

The combination of LDA and dynamical mean-field 
theory (DMFT)^^ has correctly described the Ml 
phase. However, a parameter free theory to describe cor- 
rectly both the metallic and the insulating Ml phase is 
still absent. 

The self-consistent COHSEX approximation of the 
GW-method 2 ^ gives a good description of quasi-particle 
states and has been applied to VO2 rather recently 2 ^ with 
interesting results and conclusions. Furthermore, several 
authors have used hybrid functionals, a rather ad hoc 
procedure, and this approach seems to be a good com- 
promise in many systems 2 ^ including V02>^ There is no 
good theoretical justification for the success of this pro- 
cedure. 

Optical absorption experiments create an interacting 
electron-hole pair, the exciton. Good agreement between 
theory and experiment can only be achieved by taking 
into account the exciton, especially if the system is a 
semiconductor or an insulator. Small-gap semiconduc- 
tors and metals screen excitons. For accurate absorp- 
tion spectra where such excitonic effects are important 
one has to solve the Bcthc-Salpctcr equation (BSE)2S 
which uses the intuitive quasiparticle picture. The in- 
trinsic two-particle nature of the BSE makes the calcula- 
tions very cumbersome, since a four-point equation (due 
to the propagation of two particles) has to be solved. 

In the present paper we go one step further in the 
characterization of the Ml phase of VO2 crytal. The hy- 
brid functional approach provides a better starting point, 
with the correct qualitative features, such as the correct 
symmetry of the ground state and a finite gap for a sys- 
tem in the insulating state; therefore, one expects that if 
we carry out many-body perturbation theory on top of a 
hybrid functional approach, such as HSE06 ) 30 ' 31 and we 
include the perturbative corrections, i.e., HSE06+GW 
and HSE06+GW+BSE, we will find a good convergence 
of the perturbative series. 

In the following we describe our approach in the calcu- 
lation of opto-electronic excitation of VO2 in Section HT1 
Here, we describe the applied density functionals and 
convergence tests on different parameters of the calcu- 
lations. In the next sections we present our results and 
compare them to known experimental data. First, we 
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show the calculated electronic structure with different 
methods and compare them to photo-emission data in 
Section Hill Here, we particularly analyze the phase tran- 
sition from metallic rutile structure to the Ml Mott- 
semiconductor structure of VO2. In addition, we validate 
a hybrid density functional with a parameter-free many- 
body perturbation theory method. Next, we describe 
the results on the optical properties of the Ml phase of 
VO2 in Section IIV1 again, and compare them with the 
recorded absorption spectrum. Finally, we summarize 
and conclude our results in Section IVl 



II. THE APPROACH 

We carried out DFT calculations on the Ml phase of 
VO2 as implemented within the VASP packaged We 
used small core projectors for vanadium ions, so we ex- 
plicitly included 3s3p electrons as valence. The valence 
electron states were expressed as linear combination of 
plane waves. We found that the plane wave cutoff of 
400 eV provided convergent single particle levels. We 
used the experimental geometry in all the calculations 
as shown in Fig. [TJ We attempted to optimize the struc- 
ture using the hybrid functional; the converged geometry 
disagrees with the experimental geometry rather signifi- 
cantly (the largest disagreement of the lattice constants 
was 3%), giving rise to a larger gap by a factor of 2 than 
that obtained by the same hybrid functional using the ex- 
perimental geometry. We applied various numbers of k- 
points for the unit cell depending on the calculated prop- 
erties. We found that the required size of the Monkhorst- 
Pack^ fc-point set depends strongly on the existence of 
a Mott gap. Namely, convergent charge density could be 
achieved already with a 5 x 5 x 5 Monkhorst-Pack /c-point 
set when there is a Mott-gap, while an 18 x 18 x 18 fc-point 
set was required when there was no Mott-gap. We show 
below that BSE calculations needed special treatment. 



A. The general methodology 

As a necessary step, first, we carry out standard LDA 
and GGA-type PBE^i calculations for the ground state 
density of states (DOS), the band structure, and the opti- 
cal response (the real and imaginary part of the dielectric 
function) of VO2. Then, we carry out hybrid functional 
calculations, such as HSE06 30 ' 31 which involve using a 
much more computationally demanding computational 
scheme. The reason we begin from such calculations is 
that the standard LDA or GGA calculations give a sin- 
gle particle spectrum with significant density of states at 
the Fermi level (sec Fig. and no gap. This is a qual- 
itatively different state from the correct ground state, 
and, thus, not suitable as a starting place for a pertur- 
bative scheme. The hybrid functionals, as we find, give a 
quasiparticlc gap of the same order of magnitude as the 
observed gap, and thus, we can use them as a starting 



point in a many-body perturbation theoretical approach 
to find the leading corrections. 

The HSE functional for the exchange-correlation part 
of the energy involves a parameter a which mixes the 
contribution of the short-range parts of the Hartree- 
Fock exchange and the PBE expression for the exchange 
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where the parameter uj defines what is meant by short 
and long ranged part of the Coulomb potential, which is 
split according to 
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where the part involving the complementary error- 
function erfc(uor) = 1 — erf (cor) is the short-range and 
the part involving the error-function itself is the long- 
range part. The value of the parameter u> = 0.2<2q 1 is 
determined to give a balanced description that provides 
good accuracy and speed for both molecules and solids. 31 
The E^ F ' SR (oj) is the Hartree-Fock exchange part which 
is calculated using the short-range part of the Coulomb 
interaction. Ex BE ' SR ^ LR {uj) is the Perdew et al^ ex- 
pression for the exchange energy functional which is mod- 
ified to use the short/long range part of the Coulomb 
interaction^ 

The value of the mixing parameter a = 1/4 has been 
used in a variety of materials giving rather reasonably 
accurate values for the energy gap and other quasiparti- 
cle properties. Rationale for this functional and for the 
choice of this value for a is given in Ref. HE There 
was argued that a = 1/n with n = 4 should be the 
optimum choice for typical molecules for which fourth- 
order M0ller-Plesset perturbation (MP4) yields atomiza- 
tion energies with a small mean absolute error. The 
case of n » 4 arises when there is a nearly-degenerate 
ground-state of an unperturbed problem which corre- 
sponds to the Kohn-Sham non-interacting system. An 
ideal hybrid should be sophisticated enough to optimize 
n for each system and property^ In the present paper, 
for the case of VO2 we provide results for n ~ 4 and for 
n = 8. The former is the so-called HSEOG 3 ^ functional, 
and the latter we will refer to in the following as HSE- 
1/8. We also refer to HSE06 in the following discussion 
as HSE-1/4 functional. 

As a next step, we used the so-called Go Wo approxi- 
mation to calculate the quasiparticlc spectra, as imple- 
mented in VASPi^— This means that we have used 
in Go and Wo the Kohn-Sham eigenvalues and orbitals. 
For W we take Wo = t~ x V , where the dielectric matrix 
e^ G , (q, lj), with G and G' denoting reciprocal lattice 
vectors, is calculated in the random phase approxima- 
tion and the self-energy corrections are evaluated to first 
order in the difference between the self-energy £ and the 
Kohn-Sham potential^. We note that we also carried 
out calculations where G was updated by n iterations 
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together with the DFT wave functions following the pro- 
cedure described in Ref. [39|. We denote this method by 
G„W . 

In order to calculate the optical spectra, 

a) as a first approximation, we use the independent 
particle approximation, i.e., using the energies and wave- 
functions obtained in LDA, GGA, HSE-1/4 and HSE- 
1/8. 

b) Then, we use the GW approximation on top of the 
HSE-1/4 and HSE-1/8 approximation. This corrects the 
quasiparticlc properties, such as the density of the states. 
A small correction due to the GW relative to the starting 
point indicates that the starting level of approximation 
may be good. 

c) However, transition metal oxides contain electrons 
near the Fermi level (in our case the the electrons occu- 
pying the Vanadium d-states) which are expected to be 
strongly correlated. This means that the electron-hole 
interactions are expected to be strong which should play 
a significant role in excitonic particle/hole states. These 
states are optically excited, therefore, we should and we 
will include the role of the particle/hole interaction in the 
calculation of the dielectric function. This is done using 
the Bethe Salpeter equation (BSE)<22 



B. HSE-1/4 and HSE-1/8 

In this subsection we discuss the results for the single 
particle DOS obtained with the two hybrid functionals 
HSE-1/4 and HSE-1/8 and their GW corrections. The 
DOS as a result of the GW calculation is obtained by 
finding the average energy shift (relative to the Fermi 
energy) for each quasiparticle energy level, and shifting 
the DOS calculated within DFT by the corresponding 
amount. This method is used due to the fact that a GW 
calculation with a really fine fc-space mesh is extremely 
computationally expensive. 

As can be observed in Fig. [5J the results for the DOS 
obtained using the HSE-1/4 (solid red line) functional 
and those obtained using the HSE-1/8 functional are dif- 
ferent, especially for the Vanadium d-statcs just below 
the Fermi energy. The inclusion of the Go Wo correc- 
tions has a small effect on the results of the HSE-1/4 
calculation and a much larger effect on the results ob- 
tained starting from the HSE-1/8 hybrid. Notice, how- 
ever, that the final results, i.e., HSE-l/4+GoWo and 
HSE-I/8+G0W0 are much closer to each other, especially 
for the deeper (relative to the Fermi energy) states, than 
the corresponding results without the Go Wo corrections. 

In principle, if we were to carry out a perturbation se- 
ries calculation to all orders on top of either HSE-1/4 or 
HSE-1/8, the final results should be the same. We have 
carried out higher order GW corrections, up to G4W0, 
on top of both HSE-1/4 and HSE-1/8. As expected the 
results, obtained by iterating the GW four times on top 
of the HSE-1/8, which is illustrated in Fig. [^bottom), 
show that the corrections beyond the Go Wo are small. 
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FIG. 2. (Color online) Top:Results for the density of states 
(DOS) obtained with HSE-1/4 (red solid line) and HSE-1/8 
(blue dotted line) and with HSE-l/4+GoWo (red dot-dashed 
line) and HSE-I/8+G0W0 (blue dashed line.) Bottom: Con- 
vergence study of the HSE-1/8+GW. The solid blue line is 
the result obtained with just HSE-1/8. The green dashed line 
is obtained by adding Go Wo to the HSE-1/8 and the red dot- 
dashed line is obtained after iterating GW for the 4th time 
starting from HSE-1/8. 
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Therefore, this leads us to conclude that the remaining 
difference between the HSE-1/8+GW and HSE-1/4+GW 
is not due to insufficient convergence of the GW iter- 
ation procedure. We think that this difference is due 
to processes not captured by the GW approximation. 
As we will demonstrate here by solving the BSE equa- 
tion, the effects of the quasielectron-quasihole interaction 
are large. For example, the reason for the disagreement 
between HSE-1/8+GW and HSE-1/4+GW may be the 
fact that GW does not include the role of virtual two- 
particle/one- hole excitations. Furthermore, because the 
Go Wo corrections on top of the HSE-1/4 are smaller 
than the Go Wo corrections on top of the HSE-f/8, we 
may conclude that the HSE-1/4 provides a better start- 
ing point. 



C. Convergence of BSE Calculations 

The BSE calculation is computationally demanding 
because of the fact that one has to include the interaction 
of particlc/holc pairs into different relative momenta and 
from different bands. As a result, to make the calcula- 
tion feasible within realistic computational time scales, 
we need to limit the fc-point size. 

We have studied the convergence with respect to the 
size of fc-point set used in our calculations. In Fig. [3] (top) 
we compare our calculated imaginary part of the dielec- 
tric function ti parallel to the a-axis (e2||) as calculated 
starting from the HSE-1/4 and adding the Go Wo correc- 
tion and solving the BSE for a 5 x 5 x 5 and a 7 x 7 x 7 
fc-point set ^ by including the same number of occupied 
and unoccupied bands in both calculations. All of our re- 
sults are smoothed using a exponentially weighted mov- 
ing average approach. This approach leads to the curves 
shown in Fig. [3] (bottom). Notice that the results of these 
two calculations agree reasonably well and, thus, we have 
adopted the 5x5x5 fc-point set in our calculations. 

In Fig.2]we demonstrate that using 26 occupied bands 
in our calculation maybe sufficient to achieve a satis- 
factory level of accuracy for e 2 ||- In an the calcula- 
tions presented in Fig. HI we have used 26 unoccupied 
bands (NV = 26) and we varied the number of occupied 
bands, namely, we used NO — 12, 15, 20, 26, 36. Notice 
that the results for NO = 26 have virtually converged, 
namely, these results are very close to those obtained for 
NO = 36. We have also studied the same convergence 
with respect to the number of the unoccupied bands by 
keeping the number of occupied bands fixed and found 
that when NV = 26, the results have converged. There- 
fore, all the results presented in this paper were obtained 
with NO = NV = 26. 



- — 5x5x5-N016NV8 
7x7x7 N016NV8 




— 5x5x5 


N016NV8 


-smooth 




— 7x7x7 


N016NV8- 


smooth 






Energy (eV) 



FIG. 3. (Color online) (Top) The parallel to the a-axis imagi- 
nary part of the dielectric function as calculated for a 5 x 5 x 5 
and a 7 x 7 x 7 fc-point set, by including the same number 
of occupied and unoccupied bands (18 occupied and 8 unoc- 
cupied). (Bottom) The data shown in the top part of this 
figure are smoothed using an exponentially weighted moving 
average technique. See Fig. Q] for the convention of parallel 
direction. 
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FIG. 4. (Color online) Convergence of e 2 || with number of FIG. 5. The photoemission data taken from Ref. \M is corn- 
occupied bands using a fixed number of unoccupied bands Pared with the DOS of the occupied states as obtained from 
for the case of the full BSE calculation. See Fig. [1] for the LDA, PBE calculations, 
convention of parallel direction. 



III. DENSITY OF STATES 

A. Comparison with photo-emission and 
band-structure 

In this section we present and compare the results ob- 
tained with the various levels of approximation, namely, 
LDA, PBE, HSE, and HSE+GW. In addition, we will 
compare with the photo-emission experiments. 

In Figs. El [6] we compare the density of states (DOS) as 
calculated using LDA, PBE, HSE-1/4 and HSE-1/8 with 
photoemission data4i Notice that while the LDA and 
PBE calculations give a reasonably good account for the 
occupied states which are in the range of 2 eV to 8 eV be- 
low the Fermi level, both fail to describe the states which 
lie just below the Fermi level. On the contrary, it appears 
that the HSE-1/4 calculation gives a reasonable account 
of the density of states in this narrow region. This is a 
very interesting observation, which provides hope that a 
perturbative scheme involving the GW approach and the 
BSE starting from the HSE-1/4 wave functions might be 
a good idea. 



B. Band Structure 

This hope is also justified by comparing the results of 
the band structure, obtained with HSE-1/4 with those 
obtained by GW based calculations. In Fig. the re- 
sults of Go Wo calculation on top of HSE-1/4 (solid blue 
lines) are shown as red squares. In addition, we show 
the results of G3W0 calculation with black circles. The 
overall energy scale has been changed by a constant as 



explained in the figure caption, in order for the Fermi 
levels to coincide. The overall small correction produced 
by GW indicates that the perturbative corrections on top 
of HSE-1/4 are small and, thus, we are hopeful that the 
perturbative series converges fast. 



C. Partial density of states 

In this subsection we discuss the changes that occur 
due to the phase transition from the metallic rutile phase 
to the insulating low temperature monoclinic Ml phase. 

In Fig. [8] we present the calculated partial density of 
states projected around VI and 01 (sec Fig [TJ) for the 
rutile metallic phase (top) and the insulating Ml phase 
(bottom) near the Fermi level. If we choose different 
atoms, clearly, the contribution to the partial density of 
states of different orbitals will be changed according to 
the relative position of the other atom. Notice that the 
states with the largest weight crossing the Fermi energy 
in the metallic phase are of the ti g symmetry. With 
our choice of axes (the x-axis is along the rutile c-axis, 
the y-axis is along a Vanadium- Oxide direction in one of 
the octahedra and the z-axis is along another Vanadium- 
Oxide direction in a different octahedron), the orbitals 
with t2 g symmetry are d x 2_ y 2 1 d yz and d xy , while the d z 2 
and d xz orbitals have e g symmetry. 

In the insulating phase, the states contributing the 
most in the total DOS just above the gap are mainly 
of the same ti g symmetry. Just below the gap, however, 
only the d x 2_ y 2 orbital gives most of the contribution to 
the density of states. Namely, the dimerization of the 
V-V bond splits these ti g states in opposite directions 
around the Fermi energy. The most affected orbital is 
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FIG. 6. The photoemission data taken from Ref. [4l| is com- 
pared with the DOS of the occupied states as obtained from 
HSE-1/4 and HSE-l/4+G W (top) and HSE-1/8 and HSE- 
1/8+GoWo calculations. Notice that the band gap obtained 
is about 1.0 eV with HSE-1/4 and about 0.3 eV with HSE- 
1/8. 




FIG. 7. (Color online) The band structure obtained with 
HSE-1/4 (Solid lines) is compared to that obtained by HSE- 
1/4+GoWo (red squares) and HSE-l/4+G 3 W (black cir- 
cles). The results of the GW calculations have been shifted 
downward by significant constant amounts, namely, the 
Go Wo calculation by 0.98 eV and the G3W0 calculation by 
0.60 eV to make their Fermi energies (5.50 eV and 5.18 eV 
respectively) the same with the Fermi energy of the HSE-1/4 
calculation (4.58 eV). Except for this overall constant, the 
bands from these calculations agree reasonably well with the 
HSE-1/4 calculation. 



the bonding combination of the d x 2_ y 2 orbital, one lobe 
of which is directed along the V-V dimerization direc- 
tion. Its energy is lowered relative to the Fermi level due 
to the dimerization. There is only one electron per V 
atom occupying these d ti g states which are near and be- 
low the Fermi energy. By entering the dimerized phase, 
where the unit cell doubles as shown in Fig. [TJ each of 
the two pairs of d x 2_ y 2 orbitals (one from each of the 
four V atoms in the unit cell) form one bonding and one 
antibonding linear combinations which lead to two bond- 
ing and two antibonding bands. The energy of these two 
bonding bands is separated by a gap from the two an- 
tibonding bands and from the bands formed by all the 
other tig states. Thus, the four electrons which corre- 
spond to the two pairs of vanadium atoms which partic- 
ipate in the formation of the two bonding states, occupy 
and fill up all the states in the two bonding bands. These 
correspond to the two bands just below the Fermi level 
shown in Fig. [7] This gives rise to an insulating state 
with the rest of the t^g states (along with the antibond- 
ing d x 2_ y 2 combinations) contributing to the bands just 
above the gap. 

By inspecting Fig. [71 it becomes clear that these bands 
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IV. DIELECTRIC FUNCTIONS 
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FIG. 8. The partial density of states near the Fermi level 
projected around VI for a V state and around Ol for a O 
state (see FigUJ. Top: The rutile (metallic) phase. Bottom: 
The insulating (Ml) low temperature phase. In both cases 
we have kept only the orbitals which contain more than 0.1 
electrons (out of the 4 electrons total) below the Fermi level. 



arc more or less flat (a sign of one-dimensionality due 
to the directional character of the dimerization) except 
along two directions. The direction r — > B and the di- 
rection A — > E. The former direction is along the dimer- 
ization direction (the a-axis in the Ml phase, see Fig.Q]). 
The latter direction is along the b-axis of the Ml phase 
(Fig. [T]). The former may be due to direct overlap of the 
d x 2_ y 2 orbitals along the a-axis, while the latter may be 
due to the hybridization of the d x 2_ y 2 with the oxygen 
p y . We find that the VI and V2 atoms connected to Ol 
have approximately the same distance from Ol and the 
line that begins from Ol and bisects the VI — V2 line- 
segment is our y-axis. As a result, we should expect to 
find the overlap between the p z (and the p x ) orbital of 
the Ol atom and the sum of the d x 2_ y 2 orbitals of the 
VI, and V2 atoms to should be small. This is what we 
see in Fig. [FJ where the most significant contribution to 
the PDOS just below the gap in the Ml phase is due 
to V d x 2_ y 2 and the O p y . This is not the case for the 
Ol' atom; in this case the VI and V2 atoms connected to 
Ol' are not equidistant from Ol' due to the dimerization. 
As a result, the PDOS for atom Ol' shows contributions 
from p x and p y orbitals. 



Here wc present our calculated imaginary part of the 
dielectric function for momentum q = 0, as a function 
of frequency (energy) for light polarization parallel and 
perpendicular to the c-axis of the crystal. This physical 
quantity is directly accessible by optical studies 42 

In order to calculate the dielectric function, which for 
zero momentum transfer is directly proportional to the 
optical response, we will first work within the indepen- 
dent particle approximation. In this case the dielectric 
function is calculated using the energy eigenvalues and 
eigenvectors as obtained by LDA or PBE or HSE and 
the non-interacting "bubble-diagram" for the response 
function. 

In Fig. [9J the q = dielectric function as calculated 
using LDA and PBE is compared to the experimental 
results. Notice that while the higher energy peak is in 
reasonable agreement with the experiment, this calcula- 
tion misses entirely the lower energy peak near ~ 1 eV 
and the overall distribution of strength of the response 
function. 

The electron-electron interaction, however, affects the 
excited electron-hole pair, and it can create excitonic 
effects. In the simplest extension, we can take into 
account these interaction effects by solving the Bethe- 
Salpeter equation^ In doing so, the two-particle nature 
of the the BSE equation, makes the calculation cumber- 
some, because the four point Green's function has to be 
solved self-consistently. This calculation requires very 
large memory, especially when one tries to include many 
occupied and unoccupied bands and a larger fc-point set, 
and a lot of CPU time. 

Wc found that the first allowed transition appears at 
about - 0.4 eV and - 0.2 eV for HSE-1/4+GW+BSE 
and HSE-1/8+GW+BSE methods, respectively, which 
can be tentatively compared to the calculated single par- 
ticle gaps of 1.0 eV and 0.4 eV. This is a very large re- 
duction (particularly for the HSE- 1/4 functional) which 
implies strong electron-hole interactions. We emphasize 
here that very deep occupied states could seriously in- 
fluence this value indicating the strong correlation be- 
tween the holes and electrons. It may be suspected 
that this extremely large electron-hole interaction may 
be not well accounted for by the perturbative BSE ap- 
proximation in this case. In the top of Fig. [TU] the re- 
sults for the imaginary part of the dielectric function 
e(q = 0,lu) for polarization parallel to the crystalline a- 
axis are presented for the two cases where (a) the HSE- 
1 /4 and HSE-1/4+GW+BSE and (b) HSE-1/8 and HSE- 
1 /8+GW-r-BSE are implemented and they are compared 
with the experimental data42 In the bottom of Fig. QJJ 
we present the results of the same calculations for the 
case of polarization perpendicular to the a-axis. Notice 
that the calculations based on the HSE- 1/4 functional 
show two main peaks at approximately 2 eV and 4 eV 
which correspond to the experimental peaks at approxi- 
mately 1 eV and 3 eV. Therefore, if the results of the full 
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FIG. 9. (Color online) The experimental data taken from Ref. [4^ for the imaginary part of the dielectric function at zero 
momentum as a function of frequency for the cases of polarization parallel (left) and perpendicular (right) to the a-axis is 
compared with the results of LDA and PBE. See Fig. [1] for the convention of parallel direction. 



HSE-l/4+GoWo+BSE arc shifted by a uniform amount 
of about 1 cV at lower energy they would agree with 
the experimental position of the peaks. The main rea- 
son for this disagreement arises from the HSE-1/4 cal- 
culation which overestimates the position of these peaks 
by approximately the same amount. Notice, however, 
that the BSE calculation brings the intensity of these 
peaks down to be comparable with the experimental in- 
tensity. Now, the results (right parts of Fig. [T0]) obtained 
by reducing the mixing parameter a in these calculations 
from its value of a = 1/4 in the HSE-1/4 functional to 
a = 1/8, show a better agreement with the experimen- 
tal results. This may be a fortuitous coincidence since 
HSE-1/4 provides more accurate quasi-particle energies 
as confirmed by GW calculations. In HSE-1/8 functional 
the wavefunctions are less localized than by HSE-1/4 
functional which enter BSE equations. That might com- 
pensate the neglect of the role of the electron-phonon 
interaction in the imaginary part of dielectric function. 
Furthermore, since we found that the electron-hole inter- 
action has a large effect on the dielectric function e(uj), 
one would need to include self-energy corrections due to 
two-particle+one-hole intermediate states. 



A. Optical Transitions 

It is important to note that the nature of the opti- 
cal transitions agrees well with previous analysis of the 
behaviour of this materia l 21 ' 41 , as the exact nature of 
these transitions has been in question for some time. By 
comparing the optical spectra (Fig. HU|) to the expanded 
PDOS in FigfTTl we can see quite clearly the transitions 
occurring at each peak in the spectra. It seems safe to 



do so using the projected density of states within the 
HSE-1/4 picture, as the BSE calculation did not largely 
shift the position of the peaks in the optical spectra, but 
only modulated the intensity of the quasiparticle ener- 
gies. The first peak in the optical spectra is entirely 
due to the excitation of the filled ti g states to the un- 
filled anti-bonding ti g states. No other transitions at 
this energy range are possible. At photon energies of 
~ 3eV, the transition from filled O 2p states to the same 
anti-bonding t^g states becomes available. The transi- 
tion from the filled valence band to the higher-energy 
V 3d(e g ) anti-bonding states becomes available at about 
4eV, adding to the height of the 2nd peak. The third 
peak at around 8eV arises from the transition between 
O 2p states and the anti-bonding V e g states. 

V. SUMMARY AND CONCLUSIONS 

We computed the ground state and optical properties 
of the Ml phase of the VO2 single crystal by means 
of density functional theory and beyond. We found 
that the HSE-1/4, i.e., the HSE06 functional provides 
a better description to the single particle levels as in- 
ferred from quasi-particle correction within GW approx- 
imation. However, the calculated gap is about 1.0 cV 
which is about 0.6 eV larger than gap inferred by room- 
temperature photocmission measurements. We suspect 
that the neglected vertex corrections in the quasi-particle 
correction and/or the electron-phonon interaction are re- 
sponsible for this discrepancy. 

In agreement with prior work, we find that the tran- 
sition to the Ml phase causes a Peierls-like distortion 
which, out of the three bands of t2 g symmetry, it affects 
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FIG. 10. (Color online) The experimental data taken from Ref. [4^ for the imaginary part of the dielectric function for 
polarization parallel (left column) and perpendicular (right column) to the a-axis at zero momentum as a function of fre- 
quency, is compared with the results of HSE-1/4 and HSE-1/4+GW+BSE (top row) and with the results of the HSE1/8 and 
HSE1/8+GW+BSE (bottom row). See Fig. [T]for the convention of parallel direction. 
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FIG. 11. (Color online) The PDOS projected as in Fig. [8] 
expanded to show the nature of transitions observed in the 
optical spectra. As before, states with small contribution are 
removed from the plot for clarity. 

the band which is mainly of V d x 2_ y 2 character. The four 
V atoms in the doubled unit cell after dimerization be- 
come two pairs and their corresponding four d x 2_ y 2 states 
lead to two bonding and two antibonding combinations 
separated by the Peierls gap. There are four electrons 
available for these bands and after the dimerization, the 
two bonding bands just below the Fermi energy become 
fully occupied, which leads to a Peierls insulator (for de- 
tails see our Sec. IIII Cj) . 

A major new finding in our paper is that the electron- 
hole interaction in these oxide materials is very strong 
which leads to significant reduction of the overall re- 
sponse function relative to the starting hybrid functional 
results. We believe that this may be due to the pro- 
nounced localized nature of the d x 2_ y 2 orbital and the 
strong electron correlation arising from the charge lo- 
calization. Occupied electron states with energy 10 eV 
below the Fermi-level have substantial effect on the cal- 
culated first absorption peak (~0.4 eV). This gives us 
hope that the screened Coulomb interaction between 
electrons and holes, which is strong in the Ml phase of 



VO2; might cause a high impact ionization rate, lead- 
ing to multi-exciton generation upon high-energy excita- 
tion. We believe that this finding could be rather general 
among TMOs where the same interaction is responsible 
for the Mott-gap in this class of materials. Thus, Mott- 
insulators could indeed be an important class of materials 
of interest for photovoltaic applications. In addition to 
opto-electronic excitations, doping and extraction of car- 
riers should be addressed in future studies in order to 
fully explore the applicability of this class of materials in 
third generation solar cells. 

On the other hand, our results imply that it is very dif- 
ficult to provide quantitatively accurate predictions for 
the optical response in the complex materials such as the 
Ml phase of VO2, even by means of very sophisticated 
methods involving many-body perturbation theory. The 
disagreement with the experiment on the location of the 
two main (low energy) peaks in the imaginary part of 
the dielectric function might be attributed to the fact 
that we have neglected vertex corrections and the con- 
tribution of two-particle-one hole states in the quasipar- 
ticle self-energy, the electron-phonon interaction. Fur- 
thermore, it might be significant the contribution to the 
response function from three-particle three-hole states in 
this strongly correlated material which is obviously not 
captured by BSE. Nevertheless, the HSE06+G W +BSE 
results on the optical response in the energy region im- 
portant for photo-voltaic applications, show a very large 
improvement relative to that obtained by independent 
particle approximations where the intensity of the re- 
sponse function by HSEO6+G0W0+BSE is in good agree- 
ment with the experiments. This finding indicates that 
HSEO6+G0W0+BSE method can be applied to semi- 
quantitatively predict optical excitations in transition 
metal oxides. 



VI. ACKNOWLEDGMENTS 

This work was supported in part by the U.S. Na- 
tional High Magnetic Field Laboratory which is partially 
funded by the U.S. National Science Foundation. 



* agali@eik.bmc.hu 

1 J. Mannhart and D. G. Schlom, Science 327, 1607 (2010) 

2 C. Cen, S. Thiel, G. Hammerl, C. W. Schneider, K. E. An- 
dersen, C. S. Hellberg, J. Mannhart, and J. Levy, Nature 
Materials 7, 298 (2008) 

' 5 A. Gozar, G. Logvenov, L. F. Kourkoutis, A. T. Bollinger, 
L. A. Giannuzzi, D. A. Muller, and I. Bozovic, Nature Ma- 
terials 455, 782 (2008) 

4 S. Okamoto and A. J. Millis, Nature 428, 630 (2004), (and 
references therein) 

5 A. Bhattacharya, S. J. May, S. G. E. te Velthuis, M. Waru- 
sawithana, X. Zhai, B. Jiang, J.-M. Zuo, M. R. Fitzsim- 



mons, S. D. Bader, and J. N. Eckstein, Phys. Rev. Lett. 
100, 257203(1 (2008) 

6 C. Adamo, X. Ke, P. Schiffer, A. Soukiassian, M. Waru- 
sawithana, L. Maritato, and D. G. Schlom, Appl. Phys. 
Lett. 92, 112508 (2008) 

7 A. Bhattacharya, X. Zhai, M. Warusawithana, J. N. Eck- 
stein, and S. D. Bader, Appl. Phys. Lett. 90, 222503(1 
(2007) 

8 M. P. Warusawithana, E. V. Colla, J. N. Eckstein, and 
M. B. Weissman, Phys. Rev. Lett. 90, 036802(1 (2003) 

9 e. a. M. Warusawithana, Science 324, 367 (2009) 
E. Manousakis, Phys. Rev. B 82, 125109 (2010) 



12 



11 N. F. Mott, Proceedings of the Physical Society of London 
Series A 62, 416 (1949) 

12 A. Zylbersztejn and N. F. Mott, Phys. Rev. B 11, 4383 
(1975) 

13 G. Andersson, Acta Chem. Scand. 8, 1599 (1954)10, 623 
(1956) 

14 J. B. Goodenough, Phys. Rev. 117, 1442 (1960) 

15 J. B. Goodenough, J. Solid State Chem. 3, 490 (1971) 

16 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 
(Aug 1980) 

17 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (May 
1981) 

18 R. J. O. Mossanek and M. Abbate, Journal of Physics: 
Condensed Matter 19, 346225 (2007) 

19 J. Hubbard, Proceedings of the Royal Society of London. 
Series A. Mathematical and Physical Sciences 276, 238 
(1963) 

20 G.-H. Liu, X.-Y. Deng, and R. Wen, Journal of Materials 
Science 45, 3270 (2010), ISSN 0022-2461, 10.1007/sl0853- 
010-4338-2 

21 V. Eyert, Ann. Phys. (Leipzig) 11, 650 (2002)Phys. Rev. 
Lett. 107, 016401 (2011) 

22 R. M. Wentzcovitch, W. W. Schulz, and P. B. Allen, Phys. 
Rev. Lett. 72, 3389 (1994)73, 3043 (1994) 

23 M. A. Korotin, N. A. Shorikov, and V. I. Anisimov, Phys. 
Met. Metallogr. 94, 17 (2002) 

24 S. Biermann, A. Poteryaev, A. I. Lichtenstein, and 
A. Georges, Phys. Rev. Lett. 94, 026404 (2005) 

25 J. M. Tomczak and S. Biermann, J. Phys. Condens. Matter 
19, 365206 (2007) 

26 L. Hedin, Phys. Rev. 139, A796 (Aug 1965) 

27 M. Gatti, F. Bruneval, V. Olevano, and L. Reining, Phys. 
Rev. Lett. 99, 266402 (2007) 

28 V. L. Chevrier, S. P. Ong, R. Armiento, M. K. Y. Chan, 
and G. Ceder, Phys. Rev. B 82, 075122 (2010), see Sup- 



plementary material. 

29 E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 (1951) 

30 J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 
118, 8207 (2003) 

31 A. F. Izmaylov, A. V. Krukau, O. A. Vydrov, and G. E. 
Scuseria, J. Chem. Phys 125, 224106 (2006) 

32 vienna Ab Initio Simulation Package (VASP) Version 
5.2.2.G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 
(1996)Comput. Mater. Sci. 6, 15 (1996) 

33 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 
(Jun 1976) 

34 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 
77, 3865 (1996) 

35 J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 
105, 9982 (1996) 

36 M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (Jul 
2006) 

37 M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (Jun 
2007) 

38 F. Fuchs, J. Furthmuller, F. Bechstedt, M. Shishkin, and 
G. Kresse, Phys. Rev. B 76, 115109 (Sep 2007) 

39 M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. 
99, 246403 (Dec 2007) 

40 M. S. Hyberstein and S. G. Louie, Phys. Rev. Lett. 55, 
1418 (1985)R. W. Godby, M. Schliiter, and L. J. Sham, 
56, 2415 (1986) 

41 S. Shin, S. Suga, M. Taniguchi, M. Fujisawa, H. Kan- 
zaki, A. Fujimori, H. Daimon, Y. Ueda, K. Kosuge, and 
S. Kachi, Phys. Rev. B 41, 4993 (1990) 

42 H. W. Verleur, J. A. S. Barker, and C. N. Berglund, Phys. 
Rev. 172, 788 (1968) 



